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Without Uniform-Electron-Gas Input 



Lucian A. Constantin 1 , J. M. Pitarke 2 ' 3 , J. F. Dobson 4 , A. Garcia-Lekue 1 , and John P. Perdew 5 

1 Donostia International Physics Center (DIPC), 
Manuel de Lardizabal Pasealekua, E-20018 Donostia, Basque Country 
2 CIC nanoGUNE Consolider, Mikeletegi Pasealekua 56, E-2009 Donostia, Basque Country 
z Materia Kondentsatuaren Fisika Saila, UPV/EHU, and Unidad Fisica Materiales CSIC-UPV/EHU, 

644 Posta kutxatila, E-48080 Bilbo, Basque Country 
4 Nanoscale Science and Technology Centre, Griffith University, Nathan, Queensland 4H1> Australia, 
5 Department of Physics and Quantum Theory Group, Tulane University, New Orleans, LA 70118 

(Dated: February 2, 2008) 

We resolve the long-standing controversy over the surface energy of simple metals: Density func- 
tional methods that require uniform-electron-gas input agree with each other at many levels of 
sophistication, but not with high-level correlated calculations like Fermi Hypernetted Chain and Dif- 
fusion Monte Carlo (DMC) that predict the uniform-gas correlation energy. Here we apply a very 
high-level correlated approach, the inhomogeneous Singwi-Tosi-Land-Sjolander (ISTLS) method, 
and find that the density functionals are indeed reliable (because the surface energy is "bulk-like"). 
ISTLS values are close to recently-revised DMC values. Our work also vindicates the previously- 
disputed use of uniform-gas-based nonlocal kernels in time-dependent density functional theory. 
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Density-functional theory (DFT) [l| provides ground- 
state electron densities and energies (or, in its time- 
dependent version (TDDFT) 0, excitation energies) for 
atoms, molecules, and solids. Because of its simple 
selfconsistent-field structure, DFT is used for electronic- 
structure calculations almost exclusively in condensed 
matter physics, and heavily in quantum chemistry. Exact 
in principle, the theory requires in practice approxima- 
tions for the exchange-correlation (xc) energy (or for the 
xc kernel) as a functional of the density. All commonly- 
used nonempirical approximations require input from the 
uniform electron gas, which is transferred to inhomoge- 
neous densities. The reliability of these approximations 
must be judged a posteriori, and there is a long-standing 
puzzle related to their reliability for solid surface ener- 
gies, with implications for vacancies and clusters 

The surface energy a is the energy cost per unit area 
to split a bulk solid along a plane. This is not only of 
technological importance but also a classic and highly 
sensitive test case for theories of exchange and correla- 
tion in many-electron systems. The simplest model is jel- 
lium, in which a uniform positive background of density 
n = 3 /47rr 3 = k F /Sir 3 terminates sharply at a plane and 
is neutralized by electrons that penetrate into the vac- 
uum. A local density approximation (LDA) calculation 
of jellium and simple-metal surface energies Q showed 
that the xc component a xc can be several times bigger 
than the total a, and stimulated work 0] that led to the 
development of more sophisticated functionals. 

There is now a ladder of nonempirical semilocal den- 
sity functionals, with each new rung corresponding to 
the addition of another ingredient for the energy den- 
sity. The first rung (LDA) [l[ predicts @, @| a positive 
a xc for jellium. The second rung or generalized gradi- 



ent approximation (GGA) 0] predicts values about 
3% smaller than LDA, while the third rung or meta- 
GGA @| predicts @ values about 2% lar ger than LDA. 
Ascent of the ladder brings steady improvement [6[ in 
the exactly-known Q exchange part a x . The random- 
phase approximation (RPA), which predicts inaccurate 
correlation energies for the uniform gas, predicts values 
for a xc about 6% larger than LDA [9(. The semilocal 
functionals may be corrected for long-range Coulomb ef- 
fects [6(, and the RPA may be corrected semilocally for 
short-range correlation [lOj, producing values 2 or 3% 
bigger than LDA. Use, in the framework of TDDFT, of a 
uniform-gas-based nonlocal xc kernel to correct RPA 
produces an RPA-like a xc , because RPA makes compen- 
sating errors for density fluctuations of large and inter- 
mediate wavevector 

The surface xc energies from all of the above methods 
disagree strongly with those from existing high-level cor- 
related methods: At r s — 4 (the bulk density of sodium 
metal), a xc is about 45% bigger than its LDA value in the 
Fermi Hypernetted Chain (FHNC//0) [H and Diffusion 
Monte Carlo (DMC) [13j calculations for jellium slabs. 
One conclusion from this might be that both existing 
semilocal density functionals (in the framework of DFT) 
and uniform-gas-based versions of TDDFT are not valid 
for predicting correlation energies, a very dissappointing 
outcome that would severely limit the practical useful- 
ness of DFT and TDDFT. In this Letter, we show that 
this is not the case, and in the process we resolve the 
controversy over the surface energy of simple metals. 
Refs. 
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14j already showed that a careful analysis of 
the DMC slab calculations might bring them into agree- 
ment with the density functional or RPA values, and 
also with surface energies extracted from DMC calcu- 
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TABLE I: The negative of the correlation energy per elec- 
tron (in mRyd) of a uniform electron gas in three and two 
dimensions, where r s is the radius respectively of a sphere or 
circle containing on average one electron. FHNC//0: Ref. 
QJ]. STLS 3D: Ref. STLS 2D: Ref. H|. DMC 

3D: Perdew-Wang parametrization of Ceperley- Alder Diffu- 
sion Monte Carlo (Ref. [HI). DMC 2D: Parametrization of 
Eq. (21) of Ref. Hg]. 



dim. 


r s 


FHNC//0 


STLS 


DMC 


3D 


2 


81.5 


91.4 


89.5 




3 


65.3 


74.7 


73.9 




4 


55.0 


64.0 


63.7 




5 


47.8 


56.3 


56.4 


2D 


1 




211 


219 




2 




155 


165 




4 




108 


113 




8 




66 


72 



lations for jellium spheres [l^, [Tg], suggesting that the 
surface energy puzzle had been solved. Nevertheless, one 
piece of the puzzle remained. Krotscheck and Kohn [l7| 
examined a collective RPA and used several xc kernels 
to correct for short-range effects. When they used an 
isotropic xc kernel derived from the unifor m g as, in the 
spirit of the TDDFT calculation of Ref. [TO] (see also 
Ref. [IH), they found surface energies very close to RPA, 
as in Ref. ; when they used an orbital-based Fermi hy- 
pernetted chain approximation (FHNC//0), correspond- 
ing to an anisotropic kernel constructed explicitly for the 
jellium surface, they found a large positive correction to 
the RPA surface energy. Because of this, they concluded 
that "The local-density approximation for the particle- 
hole interaction is inadequate to calculate the surface en- 
ergy of the simple metals" . 

Here we apply a very high-level correlated approach to 
calculate a, finding values that lie in the narrow range 
between meta-GGA and RPA, and much lower than the 
existing DMC or FHNC//0 slab extrapolations. We use 
an inhomogeneous orbital-based approach (ISTLS) [l9| 
that generalizes the Singwi-Tosi-Land-Sjolander (STLS) 
formalism j20|. ISTLS is, like the RPA it corrects, a fifth- 
rung density functional that employs the occupied and 
unoccupied Kohn-Sham orbitals. A comparison between 
our calculations, which do not use an isotropic xc kernel 
derived from the uniform gas, and the calculations of 
Ref. [llj leads us to the conclusion that the LDA for the 
particle-hole interaction is indeed adequate to describe 
simple metal surfaces and that existing DFT and RPA 
surface-energy calculations are reliable. 

For the homogeneous electron gas, the STLS approach 
made a remarkably accurate prediction of the correla- 
tion energy, as confirmed by later DMC calculations 
(see Table IJ). For an arbitrary inhomogeneous many- 



electron system, Dobson et al. [l9( used the linearity and 
time-invariance of a truncated Bogoliubov-Born-Green- 
Kirkwood-Yvon (BBGKY) equation [1H to propose what 
they called an inhomogeneous STLS (ISTLS) scheme, 
which can be written as a Dyson- like "screening" inte- 
gral equation for the density-response function: 



X (r,r» = x°(r,r'; W )+ / dr"Q{v, r"; w) X (r", r'; w) 



(1) 



where 



Q{t,t';u>) = - J dr"u°(r,v";i 



9(r", 



(2) 

Here, f°(r, r';u>) is a vector response function de- 
fined from the equation x°( r J r '; w ) = V r ' • i>°(r, r'w), 
X°(r,r',tL>) being the density- response function of nonin- 
teracting Kohn-Sham (KS) electrons (26|, and the equi- 
librium pair-correlation function g(r, r') is obtained from 
the fluctuation-dissipation theorem as follows [5(: 



s(r,r') = l- 



1 



7rn(r)n(r') 



du x(r, r'; iu) — 



6{r - r') 
n(r') 



(3) 

Equations (H])-© are solved selfconsistently, until a con- 
verged solution is obtained [27j |. 

By setting g(r, r') = 1 (the Hartree limit) in Eq. ([2|) 
and performing an integration by parts, the RPA is 
nicely recovered, as expected. In contrast, the inhomo- 
geneous FHNC / /0 does not recover the actual RPA but 
a "collective" RPA instead [T3|- We note, however, that 
the ISTLS scheme cannot be written using an xc ker- 
nel and does not satisfy the reciprocity constraint, i.e., 
X ISTLS (t,t';uj) / x /STL5 (r', r; w), although for jellium 
slabs this constraint is found to hold rather well. We also 
note that ISTLS is exact for all one-electron densities, for 
which g(r, r') = 0, and is exact in the high-density limit. 

For the evaluation of the ISTLS density-response func- 
tion, we extend the method described in the Appendix of 
Ref. [28| . We consider a jellium slab, assume that n(z ) 
vanishes at a distance Zq from either jellium edge [29(, 
and expand the single-particle orbitals 4>i(z) and the 
density-response function x(r, r';u>) in sine and double- 
cosine Fourier representations, respectively. Because the 
integral in Eq. ([3]) is slowly convergent and must can- 
cel out a space delta function, we use an expression 
for g(r,r') in terms of the Hartree-Fock pair-correlation 
function and the density-response functions %°(r, r';w) 
and x( r ; r '; w )- We then use the adiabatic-connection 
fluctuation-dissipation formula to obtain the xc sur- 
face energy from the following equation: 



d(q/k F )^ c 



(4) 



where q represents the magnitude of a wave vector paral- 
lel to the surface, and j* c is given by Eqs. (2) and (3) of 
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FIG. 1: Wave-vector analysis 7, of the correlation surface en- 
ergy for a jellium slab of thickness 7.21r s and r s — 2.07. Solid, 
thick dashed, thin dashed, and dotted lines represent ISTLS, 
uniform-gas based TDDFT (as reported in Ref. 11]), RPA, 
and LDA calculations, respectively, q is the magnitude of the 
2D wavevector (in the surface plane) of the density fluctua- 
tions. The area under each curve amounts to the correlation 
surface energy a c . (1 hartee/bohr 2 = 1.557 x 10 6 erg/cm 2 .) 



Ref. the density-response function Xq,\( z > z 'i w ) now 
being the 2D Fourier transform of our ISTLS density- 
response function xa(i% r'; w) °f a fictitious jellium slab 
of fixed density n(z) at coupling strength Ae 2 . 

If one replaces the interacting density-response func- 
tion Xq,)i( z i z 'i u ) entering Eq. (3) of Ref. [Il[ by the non- 
interacting density-response function then 
the exact o x is obtained, as in Ref. 0]. Here we focus our 
attention on o~ c , which for comparison we also calculate 
(i) in the LDA by replacing n c q X in Eq. (2) of Ref. [ll| 
by the uniform-gas correlation-hole density at the local 
density n(z), and (ii) within TDDFT by introducing in 
Eq. (3) of Ref. the TDDFT density-response func- 
Within TDDFT, the xc kernel 

(4) of Ref. 



tion fJ DF T{ Z ,z>;u>) 
entering Eq. 

(RPA) or the uniform-gas based isotropic xc kernel given 



ll| is taken to be either zero 



by Eqs. (6) and (7) of Ref. [11 1. 



Figure [T] shows the wave-vector analysis ^ of our 
ISTLS correlation surface energy a c (solid line), together 
with the corresponding wave-vector analysis of (i) the 
LDA correlation surface energy (dotted line), as obtained 
by using the Perdew-Wang (PW) parametrization of the 
uniform-gas correlation- hole density [30| . (ii) the RPA 
correlation surface energy (thin dashed line), and (hi) 
the isotropic xc-kernel based TDDFT correlation sur- 
face energy of Ref. 11] (thick dashed line). We observe 
that in the long-wavelength limit (q — > 0) both ISTLS 
and TDDFT calculations coincide with the RPA, which 
is exact in this limit, while the LDA fails badly [3lj ]. 
In the large-q limit, both ISTLS and TDDFT calcula- 
tions approach the LDA, as expected, while the RPA is 



TABLE II: LDA, RPA, ISTLS, TDDFT, TPSS 8], and re- 
cent DMC [33] (per Eq. |5jl) xc surface energies. Units are 
erg/cm 2 . The numerical grids we use for ISTLS are found to 
be inadequate even for RPA when r s > 3.28; nevertheless, our 
best ISTLS estimates for r s > 3.28 are found to be very close 
to the RPA. Values in parentheses represent extrapolations 
from Eq. (0. 



r s 


xc 
°LDA 


xc 

"rpa 


°"fs C TLS 


""TDDFT 


°TPSS 


"DMC 


2.00 


3357 


3467 


3417 


3466 


3380 


(3392± 50) 


2.07 


2962 


3064 


3026 


3063 


2983 


2993± 45 


2.30 


2019 


2098 


2072 


2096 


2034 


2039± 27 


2.66 


1188 


1240 


1227 


1239 


1198 


1197± 13 


3.00 


764 


801 


800 


797 


772 


768 ± 10 


3.28 


550 


579 


580 


577 


557 


551± 8 


4.00 


262 


278 


(281) 


278 


266 


(261±8) 


6.00 


53.6 


58 


(60.5) 


58 


55 


(53±...) 



wrong [32]. The important lesson that we learn from 
Fig. Q] is that two independent schemes: (i) our ISTLS 
approach, which does not use an isotropic kernel derived 
from the uniform gas, and (ii) the TDDFT approach of 
Ref. [HI], which uses a uniform-gas based isotropic xc 
kernel, yield essentially the same wave- vector analysis of 
a c . This supports the conclusion that the local-density 
approximation for the particle-hole interaction is indeed 
adequate to describe simple metal surfaces. 

To extract the surface energy of a semi-inhnite 
medium, we have considered three different values of the 
slab thickness: the threshold width at which the n = 5 
subband for the z motion is completely occupied and the 
two widths at which the n = 5 and n = 6 subbands are 
half occupied, and have followed the extrapolation pro- 
cedure of Ref. . In Table HH we show our extrapolated 
RPA, TDDFT, and ISTLS xc surface energies, for various 
values of r s . Our ISTLS calculations indicate that a per- 
sistent cancellation of short-range xc effects (beyond the 
RPA) still occurs, as in the case of the uniform-gas based 
TDDFT calculations of Ref. [11] . However, this cancel- 



lation is found not to be as complete as in Ref. [ll|, and 
yields xc surface energies that are slightly lower than in 
the RPA but still a little higher than in the LDA. Indeed, 
the difference between our ISTLS surface energies and 
their RPA counterparts is veryclose to the difference be- 
tween the conventional GGA [7J surface energies and the 
corresponding RPA-based GGA surface energies, thereby 
supporting the assumption made in Ref. [lfj that the 
short-range (beyond RPA) part of the correlation energy 
can be treated within the GGA. Our ISTLS calculations 
are also very close to the xc surface energies obtained 
by using the non-empirical TPSS meta-GGA xc energy 
functional || and a Laplacian-level metal-GGA (33[. 

The FHNC //0 approach yields a large positive correc- 
tion to the RPA surface energy. However, the FHNC//0 



4 



approach used in Ref. [l7| is in fact less accurate than 
STLS for the homogeneous 3D electron gas (see Tabled . 
STLS also does very well for the 2D electron gas (see 
Table [TJ) . These are reasons to prefer the ISTLS over 
FHNC//0 for the surface problems we are considering. 

The fixed- node DMC calculations reported by Acioli 
and Ceperley [l3| have been critiqued in Refs. 11 1, 14 1 
and [341 . Recent DMC calculations by Wood et al. 34 1 
suggest that the fixed-node approximation introduces an 
error that is slightly larger in the slab than in the bulk 
calculation and indicate that actual DMC surface ener- 
gies are larger than in the LDA but smaller than in the 
RPA, as occurs with our ISTLS calculations. 

The recent DMC calculations [34j report total surface 
energies for LDA orbitals (from which a xc is easily ex- 
tracted) and error bars for r s — 2.07, 2.30, 2.66, 3.25, 
and 3.94. To refine, interpolate, and extrapolate these 
values, we fit to them the physically-motivated form (ljj 

a xc (r s ) = A/[r 7 J 2 {\ + Bx + Cx 2 + Dx 3 )}, (5) 



where x = y/ '(1 + r s ) — 1. We choose typical values A — 



50, 000 erg/cm (correct 



limit) and D = 0.248 



(LDA fit), then vary B and C to minimize the sum of the 
squares of the fit deviation divided by the DMC error bar, 
finding B = 0.6549 and C = -0.511. Note from Table |TT] 
that LDA and TPSS both lie within the error bars of 
the recent DMC, while RPA, TDDFT, and ISTLS lie a 
little higher. The same fit has been made for ISTLS, with 
B = 0.7437 and C = -0.653. 

Finally, we note that a detailed analysis of the ori- 
gin of the xc surface energy brings us to the conclusion 
that this quantity is actually "bulk-like", arising from 
the moderately-varying-density region inside the classi- 
cal turning plane. (For r s = 2, only —3% of the total a xc 
comes from the region outside. This increases to —18% 
for r s = 4.) Inside, the reduced density gradient s falls 
in a range (0 < s < 1.9) found in the bulk of real solids, 
where gradient corrections to LDA exchange and LDA 
correlation tend to cancel. 

In summary, we have used a very high-level 
numerically-expensive correlated approach, the ISTLS 
method, to analyze the jellium surface energy into con- 
tributions from dynamical density fluctuations of various 
two-dimensional wave vectors. This analysis rules out 
the belief that the LDA for the particle-hole interaction 
might be inadequate to calculate the surface energy of 
simple metals. Furthermore, our calculations, which are 
reasonably close to uniform-gas based TDDFT calcula- 
tions [ll| and not far from the LDA, support the old 
idea that the xc surface energy should be well-described 
within LDA [Bj], and resolve the long-standing surface- 
energy controversy. 
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